﻿using System;
using System.Collections.Generic;

namespace LxCFD
{
    public partial class MGZone : IProvideZoneTopology
    {
        //这部分字段不是通过方法设置的，而是直接设置
        public MGCell[] Cells;
        public MGFace[] Faces;
        public IndexNumber InnerIndex;
        public IndexNumber TotalIndex;
        public List<IndexNumber> ThreadIndexs = new List<IndexNumber>();
        #region IProvideZoneTopology
        int[] IProvideZoneTopology.GetCellOFaces(int cellID)
        {
            return Cells[cellID].OFaces;
        }

        int[] IProvideZoneTopology.GetCellNFaces(int cellID)
        {
            return Cells[cellID].NFaces;
        }

        int IProvideZoneTopology.GetFaceOSideCell(int faceID)
        {
            return Faces[faceID].OSideCell;
        }

        int IProvideZoneTopology.GetFaceNSideCell(int faceID)
        {
            return Faces[faceID].NSideCell;
        }

        ONPair IProvideZoneTopology.GetONCellID(int faceID)
        {
            return Faces[faceID].GetONCellID();
        }

        IndexNumber IProvideZoneTopology.GetInnerIndex()
        {
            return InnerIndex;
        }

        int IProvideZoneTopology.GetThreadCount()
        {
            return ThreadIndexs.Count;
        }

        IndexNumber IProvideZoneTopology.GetThreadIndex(int threadID)
        {
            return ThreadIndexs[threadID];
        }

        public IndexNumber GetTotalIndex()
        {
            return TotalIndex;
        }
        #endregion
    }
    public partial class MGZone
    {
        public AMGInZone InAMGSolver;
        public double[] e_Common;

        public double[] e;
        public double[] ADiag;
        public double[] ALUO;
        public double[] ALUN;

        protected IAMGSolvableZone FinerZone;
        protected double[] e_FinerZone;

        #region 外部调用
        public void CreateArray(AMGInZone inAMGSolver, double[] eCommon)
        {
            LxTool.CheckCreateArray(ref e, Cells.Length);
            LxTool.CheckCreateArray(ref ADiag, Cells.Length);
            LxTool.CheckCreateArray(ref ALUO, Faces.Length);
            LxTool.CheckCreateArray(ref ALUN, InnerIndex.FaceIDs.GetCount());
            this.e_Common = eCommon;
            this.InAMGSolver = inAMGSolver;
        }

        public void SetFinerZone(IAMGSolvableZone fineZone)
        {
            this.FinerZone = fineZone;
            this.e_FinerZone = FinerZone.GetProlongArray();
        }
        public void Restrict_ArLUe(int taskIndex)
        {
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref double adiag = ref ADiag[i1];
                ref MGCell cell = ref Cells[i1];
                FinerZone.Restrict_Ar(out adiag, ref cell);
                e[i1] = cell.r / adiag;
            }
            //A_LU
            InnerIndex.FaceIDs.GetSeg(out i1, out endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                FinerZone.Restrict_ALU(out ALUO[i1], out ALUN[i1], ref Faces[i1]);
            }
            foreach (var b in InAMGSolver.e_FixValue)
            {
                set_ALUO(ThreadIndexs[b], 0, taskIndex);//理论上ALU不是0，e是0，这里直接设置ALU是0,不再初始化e
            }
            foreach (var b in InAMGSolver.e_ALUO0)
            {
                set_ALUO(ThreadIndexs[b], 0, taskIndex);
            }
            foreach (var b in InAMGSolver.e_Couple)
            {
                set_ALUO(ThreadIndexs[b], 1, taskIndex);
            }
        }
        public void Restrict_re(int taskIndex)
        {
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref MGCell cell = ref Cells[i1];
                FinerZone.Restrict_r(ref cell);
                e[i1] = cell.r / ADiag[i1];
            }
        }
        public void Prolong(int taskIndex)
        {
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                double e = this.e[i1];
                foreach (int finecell in Cells[i1].FineCells)
                {
                    e_FinerZone[finecell] += e;
                }
            }
        }
        #endregion
        protected void set_ALUO(IndexNumber indexNumber, double TargetValue, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ALUO[i1] = TargetValue;
            }
        }
        protected double cal_rMinusALU(ref MGCell coarseCell, double[] fromArray)
        {
            double rMinus = coarseCell.r;
            foreach (int faceIndex in coarseCell.OFaces)
            {
                rMinus -= ALUO[faceIndex] * fromArray[Faces[faceIndex].NSideCell];
            }
            foreach (int faceIndex in coarseCell.NFaces)
            {
                rMinus -= ALUN[faceIndex] * fromArray[Faces[faceIndex].OSideCell];
            }
            return rMinus;
        }
        protected double calALUProductValue(ref MGCell coarseCell, double[] fromArray)
        {
            double aluprotuct = 0;
            foreach (int faceIndex in coarseCell.OFaces)
            {
                aluprotuct += ALUO[faceIndex] * fromArray[Faces[faceIndex].NSideCell];
            }
            foreach (int faceIndex in coarseCell.NFaces)
            {
                aluprotuct += ALUN[faceIndex] * fromArray[Faces[faceIndex].OSideCell];
            }
            return aluprotuct;
        }
    }
    public partial class MGZone 
    {
        public void CalByJacobi_Step1(int taskIndex)
        {
            cal_x(InnerIndex, e, e_Common, taskIndex);
        }
        public void CalByJacobi_Step2(int taskIndex)
        {
            cal_x(InnerIndex, e_Common, e, taskIndex);
        }
        protected void cal_x(IndexNumber indexNumber, double[] fromArray, double[] toArray, int taskIndex)
        {
            indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                toArray[i1] = cal_rMinusALU(ref Cells[i1], fromArray) / ADiag[i1];
            }
        }
    }
    public partial class MGZone : IAMGSolvableZone
    {
        public void Restrict_Ar(out double A, ref MGCell coarseCell)
        {
            double r = 0;
            double atemp = 0;
            foreach (int cell in coarseCell.FineCells)
            {
                double adiag = ADiag[cell];
                atemp += adiag;
                r += cal_rMinusALU(ref Cells[cell], e) - adiag * e[cell];
            }
            foreach (int face in coarseCell.InnerFineFaces)
            {
                atemp += ALUO[face] + ALUN[face];
            }
            A = atemp;
            coarseCell.r = r;
        }
        public void Restrict_r(ref MGCell coarseCell)
        {
            double r = 0;
            foreach (int cell in coarseCell.FineCells)
            {
                r += cal_rMinusALU(ref Cells[cell], e) - ADiag[cell] * e[cell];
            }
            coarseCell.r = r;
        }
        public void Restrict_ALU(out double A_lu_o, out double A_lu_n, ref MGFace coarseFace)
        {
            double aluo = 0, alun = 0;
            foreach (int fineface in coarseFace.FineFaces_O2N)
            {
                aluo += ALUO[fineface];
                alun += ALUN[fineface];
            }
            foreach (int fineface in coarseFace.FineFaces_N2O)
            {
                aluo += ALUN[fineface];
                alun += ALUO[fineface];
            }
            A_lu_o = aluo;
            A_lu_n = alun;
        }
        public double[] GetProlongArray()
        {
            return e;
        }
        public void GetBounds(out List<int> e_fixValue, out List<int> e_free, out List<int> e_couple)
        {
            throw new NotImplementedException();
        }
    }
    public partial class MGZone : ISolvableZone
    {
        public void CalResidual(int taskIndex)
        {
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            double max_acphic = 0, max_Rc = 0;
            for (; i1 < endp1; i1++)
            {
                double acphic = ADiag[i1] * e[i1];
                double Rc = cal_rMinusALU(ref Cells[i1], e) - acphic;
                acphic = Math.Abs(acphic + Rc);//防止一开始acphic为0
                Rc = Math.Abs(Rc);
                if (acphic > max_acphic) max_acphic = acphic;
                if (Rc > max_Rc) max_Rc = Rc;
            }
            ThreadsMath.KeepMax(ref this.InAMGSolver.Max_Rc, max_Rc);
            ThreadsMath.KeepMax(ref this.InAMGSolver.Max_acPhic, max_acphic);
        }
        public void GetResidual_Reset(out double MaxRc, out double MaxacPhic)
        {
            throw new NotImplementedException();//不在这里实现
        }
    }
    public partial class MGZone 
    {
        public void Initialize_Step1_r0rp1rho(int taskIndex)
        {
            var r0 = InAMGSolver.BiCGSolver.r0;
            var r = InAMGSolver.BiCGSolver.r;
            var p = InAMGSolver.BiCGSolver.p;
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            double temp = 0;
            for (; i1 < endp1; i1++)
            {
                double Rc = r[i1] = r0[i1] = p[i1] = cal_rMinusALU(ref Cells[i1], e) - ADiag[i1] * e[i1];
                temp += Rc * Rc;
            }
            ThreadsMath.AddToSum(temp);
        }
        public void Cal_rhoPart1(int taskIndex)
        {
            InAMGSolver.BiCGSolver.Cal_r_Dot_r0(InnerIndex, taskIndex);
        }
        public void Cal_pArray(int taskIndex)
        {
            InAMGSolver.BiCGSolver.Cal_pArray(InnerIndex, taskIndex);
        }
        public void Cal_yArray(int taskIndex)
        {
            InAMGSolver.BiCGSolver.Cal_yArray(InnerIndex, ADiag, taskIndex);
        }
        public void Cal_vArray_alpha1(int taskIndex)
        {
            double[] v = InAMGSolver.BiCGSolver.v;
            double[] y = InAMGSolver.BiCGSolver.y;
            double[] r0 = InAMGSolver.BiCGSolver.r0;
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            double temp = 0;
            for (; i1 < endp1; i1++)
            {
                double vt = v[i1] = ADiag[i1] * y[i1] + calALUProductValue(ref Cells[i1], y);
                temp += r0[i1] * vt;
            }
            ThreadsMath.AddToSum(temp);
        }
        public void Cal_sArray(int taskIndex)
        {
            InAMGSolver.BiCGSolver.Cal_sArray(InnerIndex, taskIndex);
        }
        public void Cal_zArray(int taskIndex)
        {
            InAMGSolver.BiCGSolver.Cal_zArray(InnerIndex, ADiag, taskIndex);
        }
        public void Cal_tArray_omega1(int taskIndex)
        {
            double[] s = InAMGSolver.BiCGSolver.r;
            double[] z = InAMGSolver.BiCGSolver.z;
            double[] t = InAMGSolver.BiCGSolver.t;
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            double temp1 = 0, temp2 = 0;
            for (; i1 < endp1; i1++)
            {
                double st = s[i1];
                double tt = t[i1] = ADiag[i1] * z[i1] + calALUProductValue(ref Cells[i1], z);
                temp1 += tt * st;
                temp2 += tt * tt;
            }
            ThreadsMath.AddToSum(temp1, temp2);
        }
        public void Cal_xArray(int taskIndex)
        {
            InAMGSolver.BiCGSolver.Cal_xArray(InnerIndex, e, taskIndex);
        }
        public void Cal_rArray(int taskIndex)
        {
            InAMGSolver.BiCGSolver.Cal_rArray(InnerIndex, taskIndex);
        }
    }
}
//public void SetADLUArray(double[] A_Diag, double[] A_LU_O, double[] A_LU_N)
//{
//    this.ADiag = A_Diag;
//    this.ALUO = A_LU_O;
//    this.ALUN = A_LU_N;
//}
//public void CalByGS_Forward(int taskIndex)
//{
//    cal_x(InnerIndex, e, e, taskIndex);
//}
//public void CalByGS_Backward(int taskIndex)
//{
//    cal_x_GS_Backward(InnerIndex, taskIndex);
//}
//protected void cal_x_GS_Backward(IndexNumber indexNumber, int taskId)
//{
//    indexNumber.CellIDs.GetSeg(out int start, out int i2, taskId);
//    i2--;
//    for (; i2 >= start; i2--)
//    {
//        e[i2] = cal_rMinusALU(ref Cells[i2], e) / ADiag[i2];
//    }
//}